function [output,output2] = hm_mix_est(est_dt,sts,starting_values)

%set number of types
n_types=sts(1);

%set risk neutral or not
risk_neutral=sts(2);    %if set to zero, takes r=1 for crra utility
                        %if set to one, assumes risk neutrality
                        %if set to two, estimates r as a free parameter

%set constraints
alpha_zero=sts(3);      % sets alpha to zero
beta_zero=sts(4);       % sets beta to zero
kappa_zero=sts(5);      % sets kappa to zero
delta_zero=sts(6);      % sets delta to zero
gamma_zero=sts(7);      % sets gamma to zero

%settings for the EM-algorithm
n_steps = sts(8);       % number of starting points used
n_iterations = sts(9);  % maximal number of iterations
ll_diff_stop = sts(10);  % difference in log-likelihood for convergence

%import choices from experiment
dt=est_dt;

if risk_neutral==2
    n_params = 7;
    n_est_params = (1-alpha_zero)+(1-beta_zero)+(1-kappa_zero)+(1-delta_zero)+(1-gamma_zero)+3;
else
    n_params = 6;
    n_est_params = (1-alpha_zero)+(1-beta_zero)+(1-kappa_zero)+(1-delta_zero)+(1-gamma_zero)+2;
end

no_trials=17;
n = length(dt)/no_trials; % n is the number of subjects in the data

%sets with initial values
pm0_initial = ones(n_steps,n_types*n_params).*starting_values;
            
%constraints for the social and moral preference parameters
A_unc=zeros(n_params*n_types,n_params*n_types);
b_unc=zeros(1,n_params*n_types);
Aeq_unc=zeros(n_params*n_types,n_params*n_types);
beq_unc=zeros(1,n_params*n_types);
for type2 = 1:n_types
    if alpha_zero == 1
        Aeq_unc(n_params*type2-(n_params-1),n_params*type2-(n_params-1))=1;
        for type = 1:n_types
        if risk_neutral==2
            pm0_initial(:,n_params*type-6)=zeros(n_steps,1);   %initial values of alpha
        else
           pm0_initial(:,n_params*type-5)=zeros(n_steps,1);   %initial values of alpha
        end
        end
    end
    
    if beta_zero == 1
        Aeq_unc(n_params*type2-(n_params-2),n_params*type2-(n_params-2))=1;
        for type = 1:n_types
        if risk_neutral==2
            pm0_initial(:,n_params*type-5)=zeros(n_steps,1);   %initial values of beta
        else
           pm0_initial(:,n_params*type-4)=zeros(n_steps,1);   %initial values of beta
        end
        end
    end
    
    if kappa_zero == 1
        Aeq_unc(n_params*type2-(n_params-3),n_params*type2-(n_params-3))=1;
        for type = 1:n_types
        if risk_neutral==2
            pm0_initial(:,n_params*type-4)=zeros(n_steps,1);   %initial values of kappa
        else
           pm0_initial(:,n_params*type-3)=zeros(n_steps,1);   %initial values of kappa
        end
        end
    end
    
    if delta_zero == 1
        Aeq_unc(n_params*type2-(n_params-4),n_params*type2-(n_params-4))=1;
        for type = 1:n_types
        if risk_neutral==2
            pm0_initial(:,n_params*type-3)=zeros(n_steps,1);   %initial values of delta
        else
           pm0_initial(:,n_params*type-2)=zeros(n_steps,1);   %initial values of delta
        end
        end
    end  

    if gamma_zero == 1
        Aeq_unc(n_params*type2-(n_params-5),n_params*type2-(n_params-5))=1;
        for type = 1:n_types
        if risk_neutral==2
            pm0_initial(:,n_params*type-2)=zeros(n_steps,1);   %initial values of gamma
        else
           pm0_initial(:,n_params*type-1)=zeros(n_steps,1);   %initial values of gamma
        end
        end
    end  

end

%store results of each starting value
fit_subresults = zeros(n_steps,4);                  
pm_subresults = zeros(n_steps,n_types*n_params);    
pm0_subresults = zeros(n_steps,n_types*n_params);
pi_subresults = zeros(n_steps,n_types);
tau_subresults = zeros(n,n_types,n_steps);

setts = [n_types,n,no_trials];

% start EM algorithm
parfor s = 1:n_steps
   
    pis = zeros(n_iterations,n_types);
    pms = zeros(n_iterations,n_types*n_params);
    fits = zeros(n_iterations,4);
    taus = zeros(n,n_types,n_iterations);
    ind_lls = zeros(n,n_types,n_iterations);
    
    weighted = zeros(n,n_types);
    iter = 0;
    ll_diff = 100;

while (ll_diff>ll_diff_stop && iter < n_iterations)
    
    iter = iter+1;
    
    if iter==1        
            pm0 = pm0_initial(s,:);                         % starting values for preference parameters    
            pi0 = (1/n_types)*ones(n_types,1);              % initial values of phi (mixtures)
            pm0_subresults(s,:)=pm0;                 
            ind_ll0=lolik_mix_estep(pm0,dt,setts,risk_neutral);        
    else  
            pm0 = pms(iter-1,:);                            % starting values for preference parameters
            pi0 = pis(iter-1,:);                            % initial values of phi (mixtures)
            ind_ll0=ind_lls(:,:,iter-1);      
    end
            
    ind_l0=exp(ind_ll0);
            
            for gg = 1:n_types
                for nn = 1:n
                    weighted(nn,gg)=pi0(gg)*ind_l0(nn,gg);
                end
            end
            
            for gg2 = 1:n_types
                for nn2 = 1:n
                    taus(nn2,gg2,iter)= weighted(nn2,gg2)/(sum(weighted(nn2,:)));
                end
            end
            
    taus1=taus(:,:,iter);
    pi1=mean(taus1);
    pis(iter,:)=pi1;
    
    try
        f_unc=@(pm2)lolik_mix(pm2,dt,setts,taus1,pi1,risk_neutral);
        [pm2,fval] = fmincon(f_unc,pm0,A_unc,b_unc,Aeq_unc,beq_unc);
        il=lolik_mix_estep(pm2,dt,setts,risk_neutral);
    catch
        warning('Problem using function.');
        fval = 99999;
        pm2 = 99*ones(1,n_params*n_types);
        il = -99*ones(n,n_types);
    end
      
    pms(iter,:) = pm2;                                                              % estimated preference parameters
    fits(iter,4) = -fval;                                                           % approximated log-likelihood M-Step ("Q-value")
    ind_lls(:,:,iter) = il;
    ind_lik = exp(il);
    tot_loglik = sum(log(ind_lik*pi1'));            
    fits(iter,1) = tot_loglik;                                                      % log-likelihood
    fits(iter,3) = -trace(taus1*log(taus1'));                                       % entropy
    fits(iter,2) = -2*(tot_loglik)+(n_types*n_est_params-1)*log(n)+fits(iter,3);    % (approximated) ICL criterion (see McLachlan et al 2019) 
    
    if iter>1
        ll_diff = fits(iter,1)-fits(iter-1,1);
    end

    logl=fits(1:iter,1);
    [M,I]=max(logl);
    
    fit_subresults(s,:) = fits(I,:); 
    pm_subresults(s,:) = pms(I,:);
    pi_subresults(s,:) = pis(I,:);
    tau_subresults(:,:,s) = taus(:,:,I);
    
end

end

logl2=fit_subresults(:,1);
[M2,I2]=max(logl2);

fit_results = fit_subresults(I2,:); 
pm_results = pm_subresults(I2,:);
pm0_results = pm0_subresults(I2,:);
pi_results = pi_subresults(I2,:);
tau_results = tau_subresults(:,:,I2);

pm_matrix = reshape(pm_results,[], n_types);

output = tau_results;
output2 = pm_results;